Dynamical Clustering of Counterions on Flexible Polyelectrolytes 
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Molecular dynamics simulations are used to study the local dynamics of counterion-charged poly- 
mer association at charge densities above and below the counterion condensation threshold. Sur- 
prisingly, the counterions form weakly-interacting clusters which exhibit short range orientational 
order, and which decay slowly due to migration of ions across the diffuse double layer. The cluster 
dynamics are insensitive to an applied electric field, and qualitatively agree with the available ex- 
perimental data. The results are consistent with predictions of the classical theory only over much 
longer time scales. 
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Many biochemical reactions involve local interactions 
between small molecules (ligands) and biological macro- 
molecules such as DNA, which are polymers consistingof 
repeating ionizable groups known as polyelectrolytes 
When dissolved in a polar solvent such as water, poly- 
electrolytes ionize and counterions dissociate from the 
polymer, leaving an oppositely charged polymer back- 
bone. The highly ionized polyelectrolyte in solution at- 
tracts small mobile counterions of opposite charge, par- 
tially screening the backbone charge. The approach of 
the small molecules to a polyelectrolyte is governed to a 
large degree by electrostatic interactions, which in turn 
depend on the local molecular environment, and in par- 
ticular is sensitive to the instantaneous distribution of 
counterions around the polyelectrolyte. Thus, knowledge 
of the counterion distribution at atomic resolution is cru- 
cial for many aspects of DNA dynamics such as molecular 
assembly and inter-cellular transport. 

The difficulty in studying polyelectrolyte solutions re- 
sides in the delicate interplay of the long range elec- 
trostatic interactions, short range intcrmolccular inter- 
actions, and thermal energy, which are comparable to 
each other in magnitude. Furthermore, typical polyelec- 
trolytes are highly charged, which precludes straightfor- 
ward application of the usual Debye-Hiickel theory Q 
of electrolytic solutions. A key advance in understand- 
ing the equilibrium state of counterions associated with 
a charged polymer was the condensation theory of Man- 
ning [3| , which treated the polyelectrolyte molecule as an 
infinite rigid rod of uniform charge density —z p e/b, as an 
approximation to discrete groups of equal charge — z p e 
separated by a distance b, and considered independent 
dissolved counterions of charge z c e placed in a uniform 
bulk solvent of dielectric constant e r about the molecule. 
Considering the two-dimensional partition function in 
the plane normal to the polyelectrolyte axis, Manning 
found that the free counterion configuration is unstable 
for sufficiently strong electrostatic interaction, when the 
dimensionless "Manning parameter" £ = z p z c Ib/o > 1, 



where ls = e 2 / (e r kBT) is the "Bjerrum length" at which 
thermal energy ksT equals electrostatic potential energy. 
He then hypothesized that counterions would condense 
uniformly on the polyelectrolyte backbone so as to re- 
duce £ to unity, while the remaining (dilute) unbound 
ions in solution could be treated in the Debye-Hiickel ap- 
proximation. Since the seminal work of Manning, several 
analytical and numerical investigations Q have relaxed 
various assumptions in the original model. While the 
qualitative aspects of Manning's predictions survive, the 
phase transition at £=1 is present only for an infinite rod. 

In contrast to the average properties of the equilib- 
rium state of condensed counterions, little is known about 
the local dynamics of counterion-charged polymer as- 
sociation, which trigger and control ion-exchange reac- 
tions and affect the binding affinity of polyelectrolytes 
by shielding local electrostatic forces. In particular, most 
previous theoretical work is of a mean-field nature, and 
certainly cannot address any issues relating to temporal 
or spatial fluctuations in the distribution of the coun- 
terions along a polyelectrolyte. The electrophoresis of 
polyelectrolytes raises the further issue of the response of 
the charge distribution to an applied electric field; previ- 
ous treatments assume the counterion distribution is un- 
changed from its equilibrium state, but at least at high 
field strengths some modification will occur. 

In order to understand the fine scale behavior in the 
distribution of the counterions in the vicinity of a polyion, 
we employ molecular dynamics (MD) simulations [5[ to 
investigate, for the first time, the dynamics of counte- 
rion condensation on a mobile and flexible polyelectrolyte 
in a solution. Previous numerical studies using MD, or 
Langevin dynamics or Monte-Carlo methods have 
focused instead on the average behavior of the charge 
distribution or on polymer dynamics issues. Our com- 
putational formulation is similar to that of Chang and 
Yethiraj [|| and is based on a bead-spring model of a 
flexible polymer chain Q of Lennard-Jones monomers 
bound by FENE interactions, suspended in explicit sol- 
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vent molecules, some of which are charged to represent 
counterions and ions from added salt. The polyelec- 
trolyte chain contains N=50 monomers carrying a to- 
tal charge —Z p e distributed uniformly along the chain. 
For computational tractability, the solvent molecules are 
modeled as Lennard- Jones (LJ) spheres, with the polar- 
izability of water accounted for through its dielectric con- 
stant e r = 78 in the Coulomb interaction between charges. 
The polyion monomers and dissolved ions are likewise 
treated as LJ spheres, with an appropriate charge at 
the center, and for simplicity all are assigned the same 
mass m and core size a. All dissolved ions are mono- 
valent and the counterions from the salt and the poly- 
electrolyte are assumed to be identical. Accordingly, 
the system includes A_ = 100 coions of charge — e and 
N + = (iV_ + Z p ) counterions of charge +e so as to main- 
tain electro-neutrality. 

The simulations involve a total of 32,000 atomic units 
(monomers, co- and counter-ions and solvent) in a pe- 
riodic simulation cube of side L = 34.2a, giving a total 
number density p = 0.8c~ 3 . Simulations for one case 
with a larger box size 1.65L gives very similar results. 
From the density of water at room temperature, assum- 
ing for simplicity a cubic packing of spheres, one esti- 
mates a ~ 0.38nm. The temperature in the simulations 
is maintained at T = l.Qs/ks — 300 A' using a Nose- 
Hoover thermostat (e is the depth of the LJ potential) 
and gives a Bjerrum length Ib ~ 1.85a. Given A_ above, 
the salt concentration is approximately 0.073M, and the 
Debye length is Ijj = (4ttIb ^Pi)^ 1 ^ 2 ~ 2.7a, where Zi 
and pi are the valence and the number density of the 
simple ion species i. These parameters correspond to 
the typical experimental situation a < Ib < Id < Rg, 
where R g is the radius of gyration of the polyelectrolyte 
(found to be 5-10cr). To simulate the effects of exter- 
nal electric fields, a force which is proportional to the 
charge is added to every charged particle in the sys- 
tem. The natural MD unit of electric field strength E 
is e/ea ks 67kV/mm, and here we compare E = and 
1.0. We vary the polyelectrolyte charge Z p , such that 
the Manning parameter £ for our system (with b fs a, 
z c = 1 and z p = Z p /N) spans a range above and be- 
low the critical value 1. Lastly, the characteristic time 
scale in the simulations is r = a(mj e) 1 ! 2 , about 1 ps 
here. Typical simulations equilibrate for lOOOr and data 
is accumulated over a 3000r interval. 

To quantify the spatial distribution of the charge near 
the polyelectrolyte, we form a tube of radius r around 
the polymer chain backbone by uniting all the spher- 
ical regions of radius r centered at every monomer in 
the chain, and count the number of ions inside the tube. 
We consider a counterion to be "bound" to the polyelec- 
trolyte molecule whenever it is within a tube of radius / b ; 
the Bjerrum length is the relevant length scale because 
within this distance the electrostatic interaction domi- 
nates Brownian motion and could bind the charges to- 



gether. From the time-averaged probability distributions 
of the two nearest neighbor distances of a bound coun- 
terion to any polyelectrolyte monomer, shown in Fig. [TJ 
we see that the bound counterions prefer to sit between 
two adjacent monomers forming a triangle. Indeed, the 
typical separations are very close to each other, as well 
as to the mean monomer separation (~ la), and corre- 
spond to the three particles lying at the vertices of an 
equilateral triangle subject to modest thermal fluctua- 
tions. Note that the distributions sharpen as the poly- 
electrolyte charge increases, and that application of an 
electric field of even unit strength has little effect. 
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FIG. 1: Probability distribution of the nearest (r m i n i) and the 
next nearest (r m i n 2) distances of a bound counterion from a 
N=50 polyelectrolyte of charge Z p =20, 40 and 60. The inset 
shows the distributions when a unit electric field is applied. 

We can thus associate each bound counterion to its two 
nearest monomers in the polymer chain, which wc term a 
monomer- ion triplet. Nearby triplets are strongly corre- 
lated at short length scales because the positively charged 
counterions tend to avoid each other. To characterize 
these orientational correlations, we consider the angle 9 
between the normals to the planes formed by any two 
triplets; see Fig. EJa), and compute the time-averaged 
value of cos# as a function of the separation along the 
molecule, as defined by the difference Aj in the label j 
{j=l, 2,- • • , N), of the first monomer in the triplet. As 
seen in Fig. [2jb) , neighboring triplets are strongly anti- 
correlated out to a distance of about 5 monomers, while 
bound counterions separated by larger separations are 
independent. The correlations are weakly dependent on 
Z p , and not significantly affected by an electric field. 

The counterion-monomcr triplets exhibit dynamical 
clustering along the polyion. We define a cluster of size 
n to consist of n adjacent triplets along the chain back- 
bone with no intervening vacancies, and plot the proba- 
bility distributions P(n) for clusters of size n for differ- 
ent polyion charge Z p in Fig. [3] Provided n is not too 
small, we find an exponential distribution P(n) ~ —n, 
both with and without an electric field. Note that the 
range of charges < Z p < 60, corresponds to < £ < 2.225, 
and spans the Manning transition region. If we inter- 
pret this probability distribution as arising from a Boltz- 
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FIG. 2: (a) Geometry defining the relative orientation be- 
tween two triplets: 9 is the angle between the normals ni and 
ri2 to the planes formed by two triplets, (b) Correlation in 
the relative orientations between triplets as a function of the 
separation along the molecule, with and without electric field. 



mann factor, In P(n) ~ —En/ksT, where E n is the en- 
ergy of a cluster of size n, then the data indicates that 
E n is approximately proportional to n. In consequence 
E n i+ n 2 ~E n i+E n 2, or in another words, the energy of a 
cluster beyond a certain (small) size depends only on its 
size, independent of its environment. Thus, clusters do 
not interact with each other and are independent. Re- 
markably, the cluster energies E n are not significantly 
affected by the applied field, as can be seen from the 
weak dependence of the slopes of the linear regions of 
the P(n) curves for E = and E = 1. 
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FIG. 3: Probability of finding a monomer-ion cluster of size 
n for (a) £=0 and (b) E=l. 

The connection between independence of the clusters 
and the exponential form of P(n) can be seen in another 
way, by reference to the random sequential adsorption 
(RSA) model for polymer reactions. If ions are inde- 
pendently adsorbed with probability p and released with 
probability q per unit time at different sites, then the 
probability of forming a cluster of size n is proportional 
to \p/{j>+q)] n dH- The latter expression has the same ex- 
ponential decay with n as our simulation data, although 
this model is too simple to predict the slope correctly. 
If the ratio p/q is adjusted to give the mean number of 
counterions within the Bjerrum layer separately for each 
Z p then the slope increases systematically with Z p as in 
our simulations, but does not have the correct numerical 
value. 



Furthermore, the binding of the counterions to the 
chain is transient rather than permanent. The direct 
evidence comes from snapshots of successive configura- 
tions from the simulations, where one sees that bound 
counterions continuously move in and out of the Bjer- 
rum layer and do not remain adjacent to a particular 
monomer indefinitely. We quantify these observations 
through the charge autocorrelation function Cqq(^) of 
the total ionic charge Q{t) within the Bjerrum layer. As 
shown in Fig. |4] Cqq decays over a time close to the 
Debye time td and then fluctuates about zero, indicat- 
ing counterion decorrclation and transport through the 
Bjerrum layer. The relevant decorrelation time scale is 
the time td = 1 2 d/D for an isolated counterion to dif- 
fuse through a Debye length, where D is the diffusion 
constant of the counterions. The latter parameter was 
measured in independent simulations without a polyion 
to be D = 0.070<7 2 /t, leading to tq~ 0.1ns. Note in par- 
ticular that there is no significant difference between a 
neutral and charged polyelectrolyte molecule, and again 
an applied electric field has little effect. 
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FIG. 4: Autocorrelation functions C Q Q(t)=[{Q(t)Q(0)) - 
<Q(0)) 2 ]/[(Q(0) 2 > - <<2(0)) 2 ] for the net ion charge Q(t) in- 
side t=Ib for a polyelectrolyte of length N=50 (a) in ther- 
modynamic equilibrium and (b) translating in the presence of 
an applied field E=1.0e/ea. Each curve is an average over 5 
different realizations. The error bars shown correspond to 1 
standard deviation for the case Z p =60. 

Aside from the fluctuations apparent in Cqq, the spa- 
tial distribution of charge within the Bjerrum layer is 
constant on average, and shows an approximate agree- 
ment with Manning's predictions. Fig. [5] shows the total 
charge within a tube of radius r about polyions of charge 
Z p = 20, 40 and 60, over a time interval of 3000r^>r£> 
in the steady state. In the latter two cases, the linear 
charge densities of the polyion are above the critical value 
for Manning condensation. In the absence of an electric 
field, we see that the time-averaged total charge of the 
ions within the Bjerrum layer when £ > 1 approximately 
agrees with Manning's prediction. The original model 
was imprecise about exactly where the condensed coun- 
terions would be located; we find that the condensed 
charge in the region r < 1.5c is very close to the pre- 
diction. Moreover, the number of the coions within the 
Debye layer is found to be negligibly small as was as- 
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sumcd in the theory. Furthermore, the slopes of the solid 
curves for E = in Fig. [5] for the cases Z p =40 and 60 
above the transition point are almost the same for r>lry- 
This fact is qualitatively consistent with Manning's treat- 
ment, where the "unbound" counterions are treated by 
the linear Debye-Huckel theory, neglecting the effects of 
the polyelectrolyte charge. Also as indicated in Fig. [3 
the effect of the applied electric field is to reduce the 
charge in the vicinity of the polyelectrolyte compared to 
the zero-field value, in contrast to Manning's theory of 
the electrophoresis of polyelectrolytes for £ > 1 [l2| and 
its generalizations [l3, 14, 15. [lij. 
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FIG. 5: Charge distribution around a polyion for Z V —2Q, 40 
and 60 with and without an applied electric field. The hor- 
izontal arrows indicate the charge of condensed counterions 
predicted by Manning's theory in each case for £>1. 

Experiments have just begun to probe counterion 
structure and dynamics in polyelectrolyte solutions at 
nanometer and nanosecond scales. Small-angle neutron 
scattering measurements of the counterion-countcrion 
structure factor [l7| indicate that the counterions are 
strongly correlated with the polymer chain so as to mimic 
the pair correlations typical for the polymer, indicat- 
ing that counterions are localized in the vicinity of the 
polyions. Similarly, small-angle neutron and X-ray mea- 
surements [HI provide evidence for counterion localiza- 
tion below the Debye layer, and inelastic X-ray scatter- 
ing measurements [19[ find phonon-likc excitations in the 
counterion distribution, indicating the presence of a cor- 
related structure. Electron paramagnetic resonance spec- 
troscopy studies pijj of counterion dynamics gave direct 
evidence for the existence of preferred sites corresponding 
to the charged groups of the polyelectrolyte, where the 
counterions are electrostatically attached directly with 
no solvent molecule in between. Such "site bound" coun- 
terions presumably correspond to the bound counterions 
inside the Bjerrum layer observed in our simulations. We 
find that such counterions enter and leave the Bjerrum 
layer with a time scale around td ~ 0.1ns, while the ex- 
periments find that the exchange rate constant for the 
the "site-bound" counterions is significantly larger than 
10 9 s _1 , implying a lifetime shorter than Ins, consistent 
with our result. Lastly, recent experimental findings re- 
lated to polyelectrolyte electrophoretic mobility 2l| ar- 



gue that an ion pairing mechanism (Bjerrum association 
|22|) resembling triplet formation may play a role in the 
process of counterion condensation. 

In summary, our simulations provide insight into 
the previously unexplored local dynamics of counterion 
condensation on a mobile and flexible polyelectrolyte, 
both above and below the condensation transition. We 
find that strong long-ranged Coulomb interactions cause 
counterions to form weakly-interacting transient clusters 
around the polyelectrolyte, which exhibit short range ori- 
cntational order. The counterion clustering is robust to 
changes in the degree of ionization as well as to the appli- 
cation of an external electric field, and the cluster decay 
results from ion diffusion in and out of the charged dou- 
ble layer about the polyelectrolyte. Both the presence of 
clusters and their finite lifetime is qualitatively consistent 
with experimental data. The classical counterion con- 
densation theory is only consistent with the local charge 
distribution at much longer time scales. Our results pro- 
vide a framework for further theoretical and experimen- 
tal studies of the dynamics of counterion-charged poly- 
mer association at nanometer and nanosecond scales, in- 
formation essential for understanding the electrophoretic 
behavior of these molecules, as well as the interaction 
and binding of ligands to biological polyelectrolytes. 
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